model {
	for (i in 1:N){
		for (t in 1:7){
		       	# Sample from trunctaed normal
		       	y[i,t] ~ dinterval(z[i,t],cut[1:2])

			# squared residuals
			e[i,t] <- z[i,t]-h[i,t]
			p2[i,t] <- pow(e[i,t],2)

			# classify new response (for post. pred. checks)
			# and compare to observed values
			znew[i,t] ~ dnorm(h[i,t],1)
			ynew[i,t] <- (step(cut[1]-znew[i,t]) + 	
				2*step(znew[i,t]-cut[1])*step(cut[2]-znew[i,t]) +
				3*step(znew[i,t]-cut[2]))-1
			Match[i,t] <- equals(ynew[i,t],y[i,t])
			for (c in 1:3){
				d[i,t,c] <- equals(y[i,t],(c-1))
				dnew[i,t,c] <- equals(ynew[i,t],(c-1))
			}			
		}

		for (c in 1:3){
			dsum[i,c] <- sum(d[i,1:7,c])
			dnewsum[i,c] <- sum(dnew[i,1:7,c])
		}
						
		P1[i] <- sum(p1[i,2:7])
		P2[i] <- sum(p2[i,2:7])
				
		
	
		# First observation equation with scaled \xi_i
		z[i,1] ~ dnorm(h[i,1],1)
		h[i,1] <- xi[i]*scale +
			delta[1]*age[i,1] + 
			delta[2]*nkids[i,1] + 
			delta[3]*hhsize[i,1] + 
			delta[4]*ownocc[i,1] + 
			delta[5]*hsequi[i,1] + 
			delta[6]*incshr[i,1] + 
			delta[7]*perminc[i,1] + 
			delta[8]*transinc[i,1] + 
			delta[9]*div[i,1] + 
			delta[10]*uempl[i,1] + 
			delta[11]*fem[i] +
			delta[12]*nw[i] +
			delta[13]*degree[i] +
			delta[14]*alvl[i] +
			delta[15]*olvl[i] +
			delta[16]*union[i,1]+
			delta[17]*zg1[i] +
			delta[18]*zg2[i] +
			delta[19]*zg3[i] +
			delta[20]*zg4[i] +
			delta[21]*lond[i]
		
		# dynamics equation
		for (t in 2:7){
			z[i,t] ~ dnorm(h[i,t],1)
			h[i,t] <- alpha + gamma*z[i,t-1] + xi[i] +
				beta[1]*age[i,t] + 
				beta[2]*nkids[i,t] + 
				beta[3]*hhsize[i,t] + 
				beta[4]*ownocc[i,t] + 
				beta[5]*hsequi[i,t] + 
				beta[6]*incshr[i,t] + 
				beta[7]*perminc[i,t] + 
				beta[8]*transinc[i,t] + 
				beta[9]*div[i,t] + 
				beta[10]*uempl[i,t] + 
				beta[11]*fem[i] +
				beta[12]*nw[i] +
				beta[13]*degree[i] +
				beta[14]*alvl[i] +
				beta[15]*olvl[i] +
				beta[16]*union[i,t]
				
				p1[i,t] <- e[i,t]*e[i,t-1]	
		}
	}	
	
	# count observed and predicted category changes
	for (c in 1:3){
		ppp.obs[c] <- mean(dsum[,c])
		ppp.pred[c] <- mean(dnewsum[,c])
	}
	
	# Priors
	alpha ~ dnorm(0,0.01)
	for (b in 1:16) {beta[b] ~ dnorm(0,0.01)}
	gamma ~ dnorm(0.5, 0.01)
	for (d in 1:21) {delta[d] ~ dnorm(0,0.01)}
	scale ~ dnorm(0,0.01)
	for (i in 1:N) {xi[i] ~ dnorm(0, sigma^(-2))}
	sigma ~ dunif(0, 10)

	cut[1] <- 0
	cut.delta ~ dexp(1)
	cut[2] <- cut[1] + cut.delta
	
	# Residual correlation
	ARerr <- sum(P1[])/sum(P2[])
	
	# correctly classified cases
	TMatch <- sum(Match[,])/(N*7)
}


